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NUMERICAL SIMULATION OF SMALL PERTURBATION TRANSONIC FLOWS 


Summary 

This report summarizes the results of our systematic study of 
small perturbation transonic flows during the past year. In our 
research, both the flow over thin airfoils and the flow over wedges 
have been investigated. Various numerical schemes have been employed 
in the study. The prime goal of the research was to determine the 
efficiency of various numerical procedures by accurately evaluating 
the wave drag, both by computing the pressure integral around the 
body and by integrating the momentum loss across the shock. Numerical 
errors involved in the computations that affect the accuracy of drag 
evaluations have been analyzed. The factors that affect numerical 
stability and the rate of convergence of the iterative schemes were 
also systematically studied. The main results of the study are: 

a) The most important factors affecting the accuracy of the 
nxjmerical solution are the truncation errors, i.e., the errors 
involved in the difference approximations to the governing 
equations. For medium size nesh computations, the errors involved 
in the first-order flux conservative scheme may result in a 20% 
disagreement in the drag evaluations, although the numerical error 
in the solution to the difference equations, i.e., iLtI , has 
been reduced to 10 ^ and the global mass conservation satisfied 

to 1%. 

-4 

b) A convergence criterion of adequate for most 

medium mesh calculations, and further reduction in /A4>/ does 

’ max 

not Improve the results significantly. 
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c) In the calculations of transonic flow over a wedge, the choice 

of the sonic point operator near the wedge surface is crucial to 

the numerical behavior. We found that the direct application of 

Murman's conservative sonic point operator to the wedge 

calculation usually resulted in an instability. The noncon- 

( 2 ) 

aervative sonic point operator that we used in the past 
stabilizes the calculations but predicts a smaller supersonic 
region than the modified conservative sonic point operator that 
we developed in the course of the present study. 

d) The shock strength is affected by the far field boundary data, 
1^, which, in turn, affects the agreement of the drag as evaluated 
by the pressure integral around the body and by momentum loss 
across the shock. Several approximate methods of evaluating <}) in 
the far field have been discussed in this study, including a 
compromise approximation that results in good agreement for the 
drag evaluations. 

e) A second-order conservative scheme that we developed can be 

applied to transonic flow calculations containing an embedded 

supersonic expansion region with shock waves, such as flow over 

parabolic-arc airfoils. However, when applied to flow calculations 

with sharp supersonic compressions, numerical oscillations appear 

even if some second-order dissipation is added. 

( 3 ) 

f) The shock fitting scheme we proposed is tested in the 

present study. We conclude that this shock fitting scheme has 
no great advantage over a flux conservative scheme if the shock 
transition is of the supersonic-subsonic type although it does 
capture the Oswatitsch-Zierep singularity more readily. However, 
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shock fitting does Improve the solution In regions of flow 

containing a supersonic-supersonic shock by eliminating spurious 

differencing across the shock. Similar conclusions are reported 

(4) 

by Hafez and Cheng. 

g) For a crude initial guess of (j) , convergence of iterative 

procedures can be improved by choosing smaller relaxation factors 

in the initial stage of calculation^^^ * and then increasing 

the relaxation factors in the later stages of calculation. The 

(4) 

accelerating iterative scheme advocated by Hafez and Cheng was 
tested for the calculation of flow over a wedge with a non-uniform 
mesh distribution. No significant improvement in convergence 
behavior has been achieved in the present study. 

In summary, we are able to calculate small perturbation transonic 
flows efficiently, including flows with sharp expansions and compressions. 
To obtain reasonable agreement in wave drag evaluations we still have to 
rely on the first-order conservative calculations that employ a highly 
refined mesh. An attractive alternative is to use a higher-order scheme, 
such as the second-order flux conservative scheme we developed. However, 
further improvements in the numerical stability of such schemes are 
necessary in order to make them applicable to a wide class of flows. 

1. Introduction 

It is inevitable that the next generation of aircraft will take 
full advantage of the developing technologies in transonic flow. A 
factor of importance for all aircraft is their cruise efficiency, which 
Is proportional to the product of the aircraft's Mach number and the 
maximum lift-to-drag ratio that may be achieved at that Mach number. 


Increases In near-sonic cruise speeds, increases in aspect ratio, or 
decreases in structural weight that can be obtained with advanced 
transonic wing designs will provide for more efficient operation of 
aircraft provided that these gains are realized without noticeable 
drag penalties. Aircraft maneuverability, as well as flutter and 
buffet boundaries, are also of major concern in the transonic regime. 

How efficiently the aircraft will operate at transonic speeds may 
depend critically on our ability to calculate transonic flows numeri- 
cally, and how efficiently these calculations can be embedded in the 
design process. It will also depend critically on our ability to test 
and modify such designs on the basis of wind tunnel measurements. 

Current analytical, experimental and numerical design efforts have 
been reviewed recently by Boerstoel^^\ Cole^^^ and Pierpcnt.^^^ 

Bailey Jameson^^^^ and Ballhaus^^^^ have reviewed recent improve- 
ments in the computation of steady and nonsteady transonic flows. Such 

(13) 

calculations were first performed for airfoils by Murman and Cole , 
and later improved by Murman 

In the last several years the ability to compute such flows has 

been utilized in design studies, with a computer used to alter the 

profile shape until given design goals are reached. Work of this 

, (14) 

nature was initiated in the hodograph plane by Garabedian et al. 

and by Boerstoel and in the physical plane by Hicks et al.^^^^ and 

Barger et al.^^^^ The ramifications of improved designs for various 

(9) 

categories of aircraft are discussed by Pierpont , and their particular 

(18) 

effect on fuel economy is more fully detailed in Povinelli et al . 

The design studies of Hicks et al. calculate the lift and drag of a 
sequence of airfoils, each perturbed slightly from the previous one. 
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seeklng changes that Improve their performance. An important goal of 

such calculations Is to reduce the wave drag of the airfoils to a 

minimum over a range of Mach numbers with the lift, maximum adverse 

pressure gradient, surface area, etc., constrained. Consequently, it 

is important that such studies use a numerical code that predicts the 

wave drag accurately. The usual relaxation procedures used are probably 

adequate for this purpose if there are no rapid expansions in the flow, 

if the airfoil is not too highly curved where the shock wave meets it, 

and if a highly refined mesh is used (see, e.g., Murman's calculation 

( 19 ) 

for the parabolic-arc' ‘^). However, the airfoils that look most 
attractive usually have quite rapid expansions in the flow and the 
number of airfoils that are examined in an optimization procedure may 
be large, precluding calculations with a highly refined mesh. If the 
calculations involve modest errors in the evaluation of the wave drag, 
an improper optimum may be selected. Current computers are fast enough 
that calculations with highly refined mesh spacings can be used to 
resolve this dilemma in two dimensions, but the optimization of swept 
wings and additional design studies on two dimensional airfoils may 
benefit from numerical procedures that insure that the wave drag is 
computed accurately and reasonably efficiently. 

This report addresses the basic difficulty en<'ountered in the 
accurate evaluation of the wave drag, and suggests appropriate remedies 
for this difficulty as well as for other sources of error in these 
computations. We illustrate this basic point by considering the follow- 
ing simple mixed-flow problem: 


(K + ^ 


XX 



eD, 
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♦y(x.O) ■ f(x) on [0,1], |(iy(x,0) - 0 elsewhere, (P) 

2 2 

and <^-^0, asx + y '*•*, 

%rhere D is the upper-half plane. 

The problem (P) is an appropriately scaled representation of 
transonic flow past a nonlifting body. Here we deal with a perturbation 
potential (}), and consequently restrict f(x) such that f(x) < 1. The 
parameter K is the usual transonic similarity parameter (M^ "1)/ 

(6(y + 1))^^^ 

The partial differential equation of the problem (P) is nonlinear 

and admits weak, i.e., discontinuous solutions. These discontinuities 

are the irrotational small perturbation approximation to the shock 

waves that arise in the real flow. Because the perturbations are small 

and the flow Mach number close to 1, the vorticity introduced by these 

shock waves is sufficiently small that it may be justifiably neglected. 

There is, however, the question of the proper embedding of these 

discontinuous solutions in the physical framework that gives rise to 

( 20 ) 

the approximation (see, e.g., Dafermos ). This is accomplished by 

noting that (P) is the analog of the conservation of mass, and that the 

solution to the difference equation that represents (P) should insure 

that mass is conserved, not only where the flow is continuous, but also 

across any discontinuities that arise Hence the equation is 

more appropriately written in the form 

7 • M = 0, where ^ = -y(K+^)^J_-'iij. (1) 

z X y — 

The auxiliary condition, 

V • ^ = 0, where q = i ^ - <t j_, (2) 

— y X 

is needed to specify the shock Jump conditions. 



We seek solutions to (1) by noting that if (1) and (2) are Judged 
to be the basic conservation laws of the problem from which (P) derives, 
then 


‘‘‘'''■'’‘'discontinuity ' ‘ 

where [ ] indicates the jump across the discontinuity in the bracketed 
entity. From (1) and (2), we find 

(K + (3) 

where ( ) means the mean of the value of ( ) evaluated on each side of 
the discontinuity. 

This representation of shock waves insures the conservation of 
mass but does not allow any vorticity to be generated and hence requires 
a loss in the axial momentum flux across the discontinuity. This 
corresponds to the wave drag that results from the shock wave. This 
momentum loss must correspond to a force on the body and, consequently, 
to the drag given by the integral of the pressure over the body. The 
mathematical form of this statement derives directly from the perturbed 
axial momentum flux given by 


7 • N = 0, where N = j. (4) 

— — '■2x3x2y-'— xy — 

Using the divergence theorem, one can convert the surface integral of 

V • N to a line integral with the contour of integration shown in figure 

1 containing the body, the discontinuity and the far-field boundary. 

Applying the relation (3) across the'- discontinuity cd , the pressure 

Integral around the body, ab, is related to the jump in along cd by 


j ab 


4> 4> dx = - 1/12 


cd ^'^x^ 


dy. 


(5) 
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Equatlon (5) can also be derived directly from the surface Integral 
of the x-momentum equation around region D, l.e.. 


2 

V • { (p + pu ) ^ + puv j } dxdy 


Introducing the small disturbance approximation to the above equation 
and keeping terms of 0(6 ), one obtains the same relation (5) between 

the pressure integral and momentum loss along the shock. 

Various numerical schemes that solve (P) with the proper jump 
relations satisfied across the discontinuity are discussed. The 
numerical difficulties encountered in the study and the methods used 
to overcome such difficulties are detailed in the following sections. 


2. Numerical Procedures 

A. Difference approximations 

The well-known type-dependent difference approximations 
developed by Murman and Cole^^^^ and modified by Murman^^^ are 
adopted in the present study. The difference approximations are 
constructed such that in the supersonic region the domain of 
dependence is properly accounted for and the difference equation 
for (P) satisfies conservation of mass. After briefly reviewing 
this first-order scheme we discuss other schemes that we have used, 
i. First-order scheme: 

Detailed discussion of the fully conservative first-order 
scheme for transonic flow calculations can be found in Murman^^^ 
and also in Jameson. .An essential point is that at every 

computational point the coefficient K r is computed first 
to determine the type of equation for (P). If the computational 
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polnt is a subsonic point, central difference approximations 
are used for all derivatives in (P); if it is a supersonic 
point, backward difference approximations are used for the 
streamwise derivative, i.e., the x-derivative in (P). At 
subsonic-supersonic and supersonic-subsonic transitions, a 
sonic point operator and the shock point operator are used, 
respectively, to insure the exact cancellation of all internal 
mass fluxes generated by the sudden change in the difference 
approximations. For the problem (P) , the sonic point operator 
is equivalent to the central difference approximation to 

- 0, and the shock point operator is represented by the 
sum of the subsonic difference approximations and the super- 
sonic difference approximations to the x-derivatives, i.e.. 


{(K + (^ ) 


« ^ K • 

XX subsonic 


+ {(K + 


0 ) 
X 


) 

XX supersonic 


yy 


0 . 


( 6 ) 


Through the use of the above difference approximations, the 
total mass flux in a given computational region is conserved. 
Notice that the difference approximation (b) is inconsistent 
vd.th equation (P) • However, across a discontinuity, the 
differential equation (P) is of little consequence, whore 
the basic conservet ic'n law is imperative. In the region of 
smooth subsonic-suDor son ic transition, the soni»: point operator 
is consistent with the governing differential equation (P). 
Applications of the above fully conservative scheme to the 
calculations of flow over an airfoil are quite successful 
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and efficient. The shock wave embedded In the flow field 
Is represented by a sharp transition covering about three 
to four grid points, and the shock strength is correctly 
represented. Typical results for flow over a parabolic-arc 
airfoil are represented in figure 2. If this is applied 
without modification to the calculations of transonic flow 
over a wedge with a sharp shoulder, an instability occurs 
that originates in the sonic region near the wedge surface. 
Because the difference approximation at supersonic points 
is unconditionally stable and is dissipative, the source 
of instability must come from the sonic point operator. We 
used Jameson’s time-like analysis for the iterative procedure, 
applied to the sonic point operator, to show that it is no.i- 
dissipative and there is no damping term in the time-like 
equation. We believe this to be the cause of instability. 

it 

Let (p represent the intermediate value in the iterative 
process, and the value of at the "n+1" iteration and 

chat of the nth it '.‘ation, then. 




2<t) . . + p 

i.J 


A 

i.j-1 


0. 


(7a) 




n+1 

i.j 



+ w 




(7b) 


Here column (x-line) relaxation is used, with w denoting the 
relaxation factor. The subscripts "i" and " j " represent the 
conventional x- and y-grld locations, respectively. If we 
Introduce the following "t ime-oor i vat ive" for the iterative 


process , i . e . , 
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( 8 ) 

then we find 

(J>* - (l-w)At/w (9) 


Substituting (8) and (9) into the sonic point equation, one 
obtains the time-like equation for the sonic point 



(1-w) 


At/w <}) 


n+1 

tyy 


0 . 


( 10 ) 


Equation (10) has no damping in t as At -► 0 and hence there 
is difficulty with rapid expansions. Thus an initial error 
that appears at the sonic point may not decay as the iteration 
proceeds. The numerical oscillation of /A^/ at the sonic point 
for wedge calculations is a typical example of such instability. 
This instability may not appear if the flow expands smoothly 
from subsonic t-^ supersonic, such as in the calculations of 
flow over a parabolic-arc airfoil. There is a s.’.mple amendment 
to eliminate such instability, but this is at the expense of 
losing the conservative property at the sonic tra..sition. A 
combination of the backward difference approximation to K + i 

X 

and the central difference approximation to provide suffi- 

cient damping and thus control the instability. Let q = 

(K + < 0 at the sonic point; the difference equation is 

given by 



( 11 ) 


where subscript "b" refers to the backward difference 


approximation and "c" to the central difference approximation. 
Again, If (8) and (9) are Introduced to (11), one has 


q 'll 


n+1 

XX 


yy 


_ ,n+l 


Ax 


2-''' w * 


= 0 . 


( 12 ) 


Equation (12) is diffusive in t if w 1 2; therefore, iterative 
procedures with sonic point operator (11) will be stable. 
Numerical tests for flow over a wedge and flows over airfoils 
confirm the above statement. 

An alternative way of stabilizing the calculations is 
to treat the flow near the sharp wedge shoulder more accurately. 
For a subsonic freestream condition, the sonic point on the 
wedge surface must appear at the shoulder, i.e., point "b" 
in figure 1. This gives 


K + 




= 0 . 


(13) 


Equation (13) can be used to determine the value of ^ at the 
point right behind the shoulder, say, point "i+1,1". 


*^i+l,l 


i}) . -KAx . 


(U) 


At other sonic points Murman's fully conservative sonic point 
operator is used. With the above implementation, the iteration 
along the column just downstream of the shoulder becomes a 
fixed boundary value iteration in <{), instead of an iteration 
with a given (j) at one end and at the other. The simple 
modification (lA) at the wedge shoulder stabilizes the 
relaxation procedures and shows a convergence behavior similar 


- 14 - 


to the non-conservative sonic point operator (11). For 
the wedge calculations, the numerical results using (11) 
and (14) are quite different; this difference will be 
discussed later. 


li. Second-order scheme 

The first-order scheme usually provides results with 
poor resolution if only a moderate number of computational 
points are employed. To obtain better resolution for the 
calculations with a given total number of computational 
points, the natural alternative is to use a higher-order 
scheme. The simplest higher-order conservative scheme is 
the second-order scheme. Following procedures analogous to 
those for the first-order flux conservative scheme, one can 
construct a second-order flux-conservative scheme for the 
supersonic region, viz.. 


5 . 2 . 


^ = (2(j). .-((). , .- 24 ). _ .+4). o .)/2Ax+-Ax4i 
^x i»J 1-1, j 1-2, j i-3,j o xxx’ 


' XX 




The corresponding sonic point operator and shock point 
operator for the x-derivatives are: 

Sonic point operator: 


■ »*1.J - '■♦l-l.J 


) . . )/2Ax + 4 Ax^4 , 

1-2, J 3 xxx 


XX 




xxx 


(15a, b) 


2 

)/Ax + Ax 4' 


(16a, b) 


Shock point operator: 


1 . 2 , 




«l+l.J-*l,j-^*l-l,j«*l-2.J -♦1-3.3)^''“' -'‘“♦xxx- 


To Insure the correct transition, the sign of the type 
dependent coefficient, K + 4» , is calculated by the difference 
approximations used in the first-order scheme. The detailed 
construction of this second-order scheme is given in 
Appendix A. Notice that the supersonic expansion region, 
dissipative and dispersive truncation errors are of the same 
order and that in supersonic compressions the truncation 
errors are non-dissipative. The dominant error is dispersive 
and serious numerical difficulties are encountered. Procedures 
that avoid these difficulties are still under investigation. 

One of the methods which has been proved effective in the 
calculations of flow with an embedded supersonic expansion 
region terminated by a shock, is to add a conservative 
dissipation term to the governing equation, i.e., 

(K+* )<^ -<C = CAx^(k4) ) , with C i 0, (18) 

X XX yy XX x 

where k is a function which has the value of 0 or 1 to insure 
that the added term is also flux conservative. The difference 
approximation of each region is given below: 

Supersonic region: ^ 


f 
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(k^ ) “ (1/Ax)(k. (J) - k ) 

^XX X 1 ^XXj^ i-1 


- (I/Ax'*) - 34 >._j^+ 34) ^_ 2 - ♦i_ 3 > 


Sonic region: “ 0, ■ 1 


- a/)'* H*i - 2*i-i + *1.2) 


Shock region: k^ = 0, k^ » 1 


■ ^*t-2 * *i-3> 


One can change the coefficient C to control the magnitude 
of the dissipation. Normally, one should choose C to be 
0(1) to insure that the scheme is still second-order 
accurate. This scheme was tested for flow over a parabolic- 
arc airfoil at different Mach numbers using the results of 
the first-order calculations for the initial guess. We 
found the above scheme to be stable. The pressure distri- 
bution along the surface is shown in figure 2a. When we 
applied the same procedure to the flow over a wedge the 
value of C had to be increased by a factor of iU to damp the 
dispersive errors. The error, /A4>/ , remains bounded and 

decreases slowly, but the pressure distribution just behind 
the shoulder still exhibits oscillations, which in turn 
affects the shock strength during the iterative process. We 
suspect that a second, weak shock in the shoulder region may 


J 
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cause oscillation. However, Che existence (or nonexistence) 
of a second shock can not be established with current numerical 
calculations. 

ill. Shock fitting scheme: 

(4) 

We, as well as others (notably Hafez and Cheng , 

( 21 ) 

Morettl ), have advocated treating the shock wave as a 
discontinuity in the computations. As discussed in reference 
3, shock fitting not only provides a more reliable solution 
in the shock region by eliminating the large truncation errors 
due to improper differencing, but also saves a considerable 
amount of computing time by avoiding the unnecessary mesh 
refinement in the shock region. In the study of transonic 
flow with an embedded shock, we have introduced a shock-fitting 
scheme in conjunction with relaxation procedures. Detailed 
implementation of shock fitting to the relaxation procedures 
has been given in reference 2 and is summarized in Appendix 
B. A convergent solution can easily be obtained, and the 
results show the expected singularity when the shock intersects 
the curved surface, even with a rather crude mesh. We have 
compared the results of shock fitting with flux conservative 
calculations and found that both schemes provide solutions 
with well-defined supersonic-subsonic shocks; shock fitting 
provides a true discontinuity shock, while the flux conserv- 
ative schemes give shock transitions covering 3-4 grid points. 
For most of the numerical studies, the shock transition that 
covers 3 or 4 grid points is considered rather successful. 


Thus, we believe that In transonic flow problems containing 
only supersonic-subsonic shocks the flux conservative scheme 
la adequate and Is easier to Implement than shock-fitting 
scheme. For flow calculations containing a supersonic- 
supersonic shock, where a flux conservative scheme provides 
poor resolution for the shock, shock fitting is bound to be 
of practical importance. 

A caveat Is in order here. We have found that if shock 
fitting is introduced with the shock in the wrong position, 
located say, by a nonconservative calculation, the values 
of may change so that the jump conditions are satisfied 
but the shock may not move to the location predicted by a 
conservative calculation (without shock fitting). This 
"nonunique" behavior is believed to be due to the fact that 
changes in (p behind the shock, required by shock fitting, 
are accommodated erroneously in the nonlinear contribution 
to far-field boundary data rather than effecting data near 
the sonic line and hence the shock position. The use of 
shock fitting in conjunction with a fast Poisson solver needs 
to be implemented to see if this resolves the difficulty. 

B. Evaluation of far field 

All the calculations shown in this report are carried out in 

a finite computational region. That is, we avoided mapping an 

infinite domain to a finite one. Thus, an approximate formula to 

( 22 ) 

evaluate ij) on the boundaries is necessary. We used Klunker's^ 
representation , 


f(0 log [(x-o^ - Ky^l d5 



( 19 ) 


where the first term represents the source contribution due to 

the body, and the second term the nonlinear contribution throughout 

the flow field. Direct evaluation of the nonlinear contribution 

requires tedious double integration at every boundary point and 

proves to be extremely time-consuming. An alternative is to pull 

2 2 

out the kernel (x-S)/{(x-^) -K(y-n) } and evaluate it at a mean 
position such that the double integral is calculated just once and 
used for all boundary points; i.e., 



X - g 

(x-O^-K(y-n)^ 


d^dn 


X - X 

av 

(x-x )^-K(y-y 
av av 




I I dWn. (20) 

*00 

A good approximation for (x^y» ^ symmetric problem is 

( -I-, 0 ), i.e., at the center point of the chord. One can also 
obtain a higher order approximation by a series expansion of the 
kernel with respect to i^n/x - i/Ky and ^ + »iCn/x + »^, i.e., 

2 2 

(x-5)/{ (x-C)"- K(y-n)^} = l/(x^-Ky^){x + ^ 

X -Ky 


2x 


2^2 
X -Ky 


2Kxy 

x^-Ky 


2 


2Kv 


2 

X 



n 


Cn + . . . } 
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In the above representation, one has to evaluate 




dCdn and 

X 


2 f f 2 

d^dn In addition to dCdn< 

X j J X 


These integrations can be obtained with little extra computational 

2 


'd^dn. We have used 


effort over that required to evaluate 

J * 

the above expansion for the calculations of far field 4>, and 

found its effect on the results is negligible if the size of the 

computational region is large enough, say, 5 by 5 chord lengths. 

But it is important that x and y be chosen rationally. This 

av av 


fact is illustrated in section 3. 


C. Convergence criteria 

The iterative procedures terminate when the numerical 
solution to the difference equations reaches a prescribed accuracy. 
There are many ways of defining numerical accuracy. The most 
common one is to set the maximum change in throughout the ccmpu- 
tatlonal field at each successive iteration to be less than a 
small quantity, e, i.e., most of the calculations, 

an e of 10~^ or 10 ^ is enough. Further reduction in c gives a 
more accurate solution to the difference equations but not necessarily 
to the differential equation unless finer meshes are used. 

Another way of defining convergence is to require that the 
residual of the difference equation at each computational point 
be less than a given small number. Normally, for given accuracy 
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9 

in • the residual of the difference equation Is at least 

max 

one order of magnitude larger than /A^/ . Thus, using the 

max 

maximum residual in the field as a convergence criterion is more 
restrictive than computing ^ to the same accuracy. In this 
report, we use /A(^/^^ < 10 ^ unless specifically mentioned. 

D. Relaxation factors and acceleration schemes 

The convergence behavior of the relaxation procedures 
depends on the choice of relaxation factors. In some cases, an 
improper guess of relaxation factors causes the scheme to diverge. 
Because the problem (P) is nonlinear, there is no simple way of 
finding optimal relaxation factors. The choice of the best 
relaxation factors usually relies on numerical experiments. As 
a general rule it has been found that it is better to choose 
smaller relaxation factors in the initial stage of the calcu- 
lations, including both the subsonic over-relaxation factor and 
the supersonic under-relaxation factor. In the later stage of 
calculations, one should increase the relaxation factors to as 
large a value as possible (within the limit of stability criteria) 
to enhance the convergence. Normally, the upper limits of the 
relaxation factors for medium mesh calculations are 1.95 for the 
subsonic region and 0.90 for the supersonic region. 

For nonuniform mesh calculations, little is known about the 
choice of relaxation factors. We have experienced a slow rate 
of convergence for a range of relaxation factors. In most cases, 
for a given relaxation factor, the rate of convergence for non- 
uniform X- and y-mesh distribution is slower than that of uniform 
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mesh distribution. Further investigations are necessary in order 
to improve the efficiency of the relaxation procedures. 

The accelerated iterative scheme proposed by Hafez and 

Cheng has also been tested on the calculations of flow over a 

wedge using non-uniform mesh distribution. We were not able to 

achieve the fast rates of convergence reported by Caughey and 
(23) 

Jameson. The extrapolated value <)) for each accelerated step 
usually increases the errors in 4> rapidly, and the relaxation step 
following the accelerated step decreases the errors to the previous 
level. Such behavior repeats in the Iterative process, which in 
turn provides little gain over a pure relaxation scheme. We believe 
such behavior is caused by the poor estimation of the dominant 
eigenvalue for the iterative matrix using a non-uniform mesh 
distribution. The difficulty of finding an accurate estimate of 
the dominant eigenvalue may restrict general application of 
accelerated iterative schemes. 

Recently, we have tried the acceleration scheme proposed by 
Jameson^^^\ i.e., a Poisson solver step followed by a relaxation 
step. The Poisson solver computes the subsonic region efficiently, 
where the relaxation step determines the numerical solution in the 
supersonic region. This improves the rate of convergence at least 
of one order of magnitude. We found the scheme works fairly well 
in the physical plane for the calculations of the small perturbation 
equation. The Poisson solver we used was obtained from NC/VR and 
solves Poisson's equation for a uniform grid. In our test 
calculations, each Poisson solver step is followed by a few 
relaxation steps, say 5 to 10 steps. Detailed discussion on such 
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a scheme will be submitted in a report on a subsequent grant. 

Typical convergence rates for the combined scheme are shown in 

-6 

figure 7. The solution converges to < 10 within 100 

max 

relaxation steps and 20 Poisson solver steps. 

3. Results and Discussions 

Most of the present numerical results were obtained using uniform 
grid calculations. In the case of flow over a sharp wedge, highly non- 
uniform mesh distributions along both x- and y-direction have also 
been used. The smallest meshes used ''.n the wedge calculations were 
Ax “ 0.01 in the nose and shoulder regions, and Ay ■ 0.01 near the 
wedge surface. The mesh Increases exponentially away from the nose, 
the shoulder, and the wedge surface. The region of computation for 
non-uniform mesh distribution is, in general, much larger than that 
of uniform mesh distribution. The typical size of the computational 
region for a uniform mesh distribtuion is 6 by 6 chord lengths, while 
for a non-uniform mesh distribution it is 25 by 15 chord lengths. In 
the airfoil calculations, it has been reduced to 3 by 3 chord length.s 
when the embedded supersonic region is relatively small. 

Figure 2a depicts the pressure distributions Cp = -2(tC 
along the parabolic-arc at M«, * 0.909. The solid line is the result 
of a first-order scheme with Ax ■ Ay ■ 0.05. The points are the 
results of a second-order scheme using the same mesh size. Because 
of the smooth supersonic expansion, there is no significant difference 
in pressure distribution away from the shock. Figure 2b shows the 
pressure distributions for the same airfoil at M® * 0.9, using both 
fine mesh (solid line) and coarse mesh first-order calculations. 
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Agaln, Che results are quite similar, except Che fine mesh calculations 
give a sharper shock. The wave drag, calculated '''om c* e results of 
each scheme for M„ - 0.909, are listed in Che following table, where 
subscript "PI" refers to the drag evaluated by pressure integral around 
the body, and "SK" by the momentum loss across the shock. 

Table 1. Comparison of wave drag for parabolic-arc at 
M, - 0.909 


Scheme 

S(PI) 

S(SK) 

1st order 

0.02875 

0.02888 

2nd order 

0.02886 

0.02918 

Murman^^^^ refin mesh 

0.0315 

0.0320 


The agreement between and within each scheme is extremely 

good for this particular case, i.e., approximately 1%. However, using 
the same mesh distribution, we computed the drag for = 0.9, and 
found and differ bv 12.5%, (0.0214 vs. 0.0245). The 

disagreement is reduced to 5% (0.0229 vs. 0.0240) by halving the x-mesh 
spacing. Thus, we conclude that the drag evaluation is rather sensitive 
to the numerical approximation to the solution of the differential 
equation, i.e., to the truncation errors of the difference approximations. 
This will be discussed further in another report. 

Figure 3 illustrates the supersonic region for flow over a wedge 


at K ■ -0.5. Curve 1 shows the sonic line and shock obtained by using 
the non-conservative sonic point operator and curve 2 that for the fully 
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conservacive sonic point operator. Because of the sharp expansion 
around the shoulder, the results are extremely sensitive to the change 
of difference approximations to the differential equation In that 
region. The non-conservative sonic point operator pushes the sonic 
line approximately 2 meshes ahead of the shoulder, whereas the modified 
conservative sonic point operator (equation 13) predicts the correct 
sonic condition at the shoulder. Inaccuracy in locating the sonic 
line may cause an lncor*ect expansion near the shoulder, which, in 
turn, results in an incorrect shock position. The pressure distributions 
along the wedge are shown in figure 4. Because conservative shock point 
operators wrre used for both calculations, the jump across the shock 
satisfies the normal shock relation even though the shock posltic.is are 
quite different. The drag values for each case are listed in Table 2. 

Table 2. Wave drags for flow over wedge at K ■ -0.5. 


Scheme 

S(pi) 

S(SK) 

Non-conservative 
sonic point operator 

0.6535 

0.5196 

Conserve 've 
scheme 

0.7627 

0.6577 

Conservative scheme 
with modified B.C. 

0.7561 

0.7577 


The flux conservative scheme gives slightly tetter agreement between 


^D(PI) ^D(SK) non-conservative sonic point calculations. 

The cause of poor agreement between and , as mentionea 

before, is due to truncation errors. In wedge calculations, we found 
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thaC Is Insensitive to the Iteration process, l.e., 

approaches Its final value much faster than reason Is 

that flow on the wedge surface that contributes to drag evaluation Is 
subsonic. The inaccuracy In the supersonic r'‘gion has little effect 
on Its pressure distribution. In contrast, is affected by the 

numerical accuracy in the supersonic region. Thus, for the wedge, 
S(PI) reliable than 

In the calculation of flow over a wedge, the truncation errors 
around the shoulder cannot be reduced unless a much refined mesh is 
used. The shock is improperly calculated because of the inaccurate 
expansion near the shoulder. One way of obtaining agreement between 
the two evaluations of wave drag is to modify the far-field boundary 
data. The boundary data for i}i can be redistributed by evaluating the 
kernel in (19) at different positions. Because is insensitive 

to iteration, one can expect that a slight modification of (19) will 
not change the pressure distribution around the wedge significantly. 
However, is sensitive to the iteration. Therefore, a redistri- 

bution of boundary data can result in a significant change in shock 
strength. We have shifted the kernel successively tu (-0.25, 0) and 
found a good agreement in and shown in Table 2. 

Interestingly enough, the global conservation of mass has not been 
affected by the modification of boundary data , i.e., the unbalanced 
flux is still within 17 .. This simply illustrates the sensitivity of 
the results to the far-field boundary data. The choice (-0.25, 0) is 
not consistent with the evaluation of the full integral (19). We 
conjecture that the large truncation errors in (4) are compensated for 
by this erroneous evaluation on the boundary. 
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To confirm that the disagreement in and is due to 
large truncation errors near the sharp shoulder, we have carried out 
the same calculations for the case of a smooth shoulder, i.e.. 


♦y ■ C, for 0 < X < 0.5, 


■ C cos {ir(x-0.5)}, for 0.5 < x < 1, 


( 22 ) 


with C » 2 it/ (tt + 2) such that the thickness of the wedge is equal to 1. 
The pressure distribution for a smooth wedge is shown in figure 5. The 
discrepancy between and is reduced to 6Z (i.e., 0.5185 

vs, 0.4845). One can also shift the kernel in (19) to (0.25, 0) instead 
of (-0.25, 0), and obtain an excellent agreement in snd 

(i.e., 0.5142 vs. 0.5106). 

The pressure distribution along a sharp wedge using highly 
non-uniform mesh calculations is shown in figure 6. Except for a 
sharper compression downstream of the shoulder, the general feature 
of the pressure distribution is quite similar to that of uniform mesh 
calculations . 

The results of a shock-fitting calculation have been included in 
the previous report and also presented in reference 2, and will not be 
repeated here. Shock fitting does provide better resolution in the 
shock region. However, the discrepancy in and in the 

wedge calculations, which is due mainly to poor resolution near the 
expansion region (the shoulder) cannot be improved by shock fitting. 

In figure 7 we showed the convergence behavior for both sharp \ edge 

and smooth wedge calculations. For a pure relaxation scheme the smoother 

flow has better convergence. Typically, it requires 500 iterations or 

more to bring the error /A*/ down to 10 

max 



Figure 2a. Pressure distribution along a parabolic-arc airfoil at 

Mo, ■ 0.909. : first-order scheme. second- 

order scheme. : with shock, fitting. See Table 1. 
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Figure 2b. Continued . Mu = 0.900. First-order scheme, : 

mesh; "x": crude mesh. See text following Table 1, 
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Flgure 3. Sonic line and shock for flow over wedge at K = -0.5 
Curve 1: non-conservative sonic point operator; 

Curve 2: "ipdified conservative scheme. 



Figure 4. Pressure distribution along wedge at K = -0.5. Curve 1: 
non-conservative sonic point operator; Curve 2: modified 

conservative scheme; "x': conservative scheme with modified 

far-field boundary condition (note change in shock strength^. 

OV'K-.lSA'. See Table 2. 

OK k'OilU 



Figure 5. Pressure distribution alon? a smooth wedge at K = -0.5. 

: conservative scheme; "x" : conservative scheme with 

modified far-field boundary condition. 
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Appendlx A ; Construction of a Second-Order Conservative Scheme 

The basic requirement for a numerical scheme to be flux 
conservative is that the difference equation be in divergence form. 
The difference approximations (15a, b) satisfy the above requirement 
when we substitute them into (K + )<j> , i.e., 

X XX 


2Ax^^^ Ax ; (.K + f 

~ lix '^x^ i-1/2 ■ i-3/2^ 


where (K + 


K + 

Ax 


It can easily be seen that for purely supersonic flow, the flux 
2 

(K + A ) at the interface cancels exactly when the difference equations 

X 

are summed with respect to i and j . This is also true for a purely 
subsonic region using central difference approximations, where 


^ i-1/2 


K + («^ 


i.j "i-l,j 


J/Ax. 


Because of the change of difference approximations at the subsonic- 
supersonic interface, special care has to be taken in forming the 
difference approximations in order to maintain the flux conservative 
property. Let us consider a strip of computational points, where "a" 
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denotes the sonic point, and "s" the shock point. For slopliclty, 
let us assume that the supersonic region is confined in a < i < s. 


M 

s 

X 


a 

s 

a 

a-1 

a 

a+l 


s-l 

s 

s+l 


Sketch A1 


The diff erence approximations for (K + 3t each point are; 


Point 


Difference Approximation for (K + ‘^x^'^xx 


a-1 




a+l 


a- 1) 


Ax 


2<i -34s .+4 T - 

(■< " ^ i:' ) > 


s-l 


- 2(ti ,- 34 > -+4 .2 2^ -3^ , : 


s+l 




s-t-2 ^s-H 
Ax 




2 

At point "a-1", the incoming flux {K + ^'^a-l~*'a-2^ ^ cancelled 

by the outgoing flux at point a-2, the flux that remains uncancelled is 
2 

{K + (4 -4 ,)/Ax} , which has to be cancelled by the difference 

a a-1 

approximations at "a". Similarly, the uncancclled flux at point 


"a+l" is 
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- (K 


Ax 



Therefore, the only choice of difference approximations for 
at point "a" that conserves mass Is 




- L -- ) - C' * > 




2Ax 


Ax 


which Is given by (16a, b). 
Similarly, at shock point "s". 


{(K+^ )<J> } 

^X XX s 


2S^ '(■' 


s+1 


Ax 


-.2 , 2(i) ,-3(J: 

4 - (K + -) 1 




2Ax 


Ax 


To insure the correct transition, the type dependent coefficient 
K + is tested at each computational point by 

qc - <K + - K + (♦i,i-»i.p/2ix 

- (K + ♦ )u ■ K + (♦,-♦, ,)/2ix, 
b X b 1 1“/ 

where the sonic point is defined by > 0, and < 0, and the shock 
point by q^ < 0, q^^ > 


0 . 


Appendix B ; Shock Fitting 


A shock fitting scheme can easily be Implemented in the relaxation 
calculations. We introduce shock fitting when a reasonably well- 
converged first-order result is obtained. The initial shock position 
Is determined by interpolating the sonic position in the i ;j>ersonic to 
subsonic transitions, and the shock points are then treated a'’ regular 
computational points. (Sketch Bl) . At the shock point, the flow 
properties ahead of the shock, i.e., 
can be determined either by 
extrapolation from the upstream 
conditions at Q and R, or by using 
the characteristic relations along 
c"*^ and c”. Both methods have been 
tested and have given satisfactory 
results. For most of the calculations, Sketch Bl 

the simple extrapolation method was used. Behind the shock, is 
extrapolated by the usual difference method using the old value at P, 
and $ is obtained from the jump relation (3). Should the shock slope 

f 

extrapolation provide a shock intersection at s as Indicated by the 
dotted line, then the value of is fixed at s' by using the value of 

^ ♦ at s. At each stage of the calculation we correct the 

X y 

position of the shock by the simple procedure 



n+1 


+ kAx (4i, 


n+1 


(B-1) 


Various values of k have been tried; values near one seem to be the 
roost satisfactory. The rationale for (B-1) Is simple: if the flow 
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ahead of the shock has increased in speed, i.e., ^ 

(B-1) allows the shock to be swept downstream and flow properties can 
be recalculated using the new shock position. 

Because the shock is treated as a discontinuity, the flow ahead 
of and behind the shock is continuous. Thus, the differenr. approxi- 
mations for the grid points adjacent to the shock point, e.g., P and Q 
can be constructed using the usual difference approximations without 
taking into consideration the exact internal flux cancellation. 
Numerical tests showed that if the difference approximations at P 
and Q are consistent with the differential equation, the internal 
fluxes generated at those points are small. 



